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© Method and device for adapUvely estimating a transfer function of an unknown system. 

• 

© In an adaptive estimation of a transfer .function by a projection algorithm, a forward linear prediction 
coefficient vector a(k) of an input signaUx(k), the sum of forward ^'posteriori prediction-error squares F(k), a 
backward linear prediction coefficient vector b(k) oi the input signal x(k) and the sum of backward a posteriori 
prediction-error:, squares B(k) arc cOmputeiJ.. Letting a step; 4\ze and a pre-fflter deriving' coefficient vector be 
r presented by u and f(k). respectively, a pre-frlter coefficient vector g(k) is calculated by a recursion formula for 
the pre-filter coefficient vector g(h). which is composed of the following first and second equations: 
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BACKGROUND OF THE INVENTION 

The present invention relates to a method for adaptively estimating, with a projection algorithm, a 
transfer function of an unknown system and its output in an acoustic canceller, active noise control or the 

5 like and an estimating device using such a method. 

In the following description, time will be represented by a discrete time k. For example, the amplitude of 
a signal x at time k will be expressed by x(k). Fig. 1 is a diagram for explaining the estimation of the 
transfer function of an unknown system. Reference numeral 11 denotes a transfer function estimation part 
and 12 the unknown system, and reference character x(k) represents an input signal to the unknown system 

10 and y(k) an output signal therefrom. The transfer function h(k)tof the unknown system is estimated using 
the input signal x(k) and the output signal y(kj. Fig. 2 is- a diagram for explaining an adaptive estimation of 
the transfer function. Reference numeral 21 denotes an estimated transfer function correcting vector 
calculation part, 22 an estimated transfer * function correction part and 23 a convolution part, these part 
constituting an estimated signal generation part 20. The transfer function h(k) of the unknown system 12 is 

75 estimated as a transfer function h(k) of an FIR filter of a tap number L which forms the convolution part 23. 
More specifically, coefficients h\(k), .... h L (k) of the FIR filter are estimated. Let it be assumed that the 
"transfer function" and the "FIR filter coefficient" will hereinafter be construed as the same. For the sake of 
brevity, the filter coefficient is represented as an estimated transfer function vector h(k) defined by the 
following equation. 

20 

h(k) = [hi(k), h 2 (k), h L (k)] T (1) 

where T represents a transpose. 
25 In Fig. 2, the input signal x(k) to the unknown system 12 is fed to the convolution part 23 and a 
calculation is performed to obtain the estimated transfer function vector h(k) that minimizes an expected 
value which is the square of an error signal e(k) available from a subtracter 24 which detects a difference 
between an output y(k) from the convolution part 23 given by the following equation (2) and the output y(k) 
from the unknown system 12. 

30 

y(k) = h(k) T x(k) (2) 
x(k) = [x(k),x(k-1), ...,x(k-L + 1)] T (3) 

35 

where y(k) is an estimated value of the output from the unknown system, which is close to the value of the 
output y(k) when the estimated transfer function h(k) is close to an unknown characteristic. 

In practice, the transfer function of an unknown system often varies with time as in the case where the 
transfer function of an acoustic path varies with movement of audiences or objects in a sound field such as 

40 a conference hall or theater. For this reason, the estimated transfer function h(k) of the unknown system 
also needs to be adaptively corrected in accordance with a change in the acoustic path of the unknown 
system. The estimated transfer function correcting vector calculation part 21 calculates an estimated 
transfer function correcting vector 5h(k) on the* basis of the error signal e(k) and the input signal x(k) to the 
unknown system 12. The estimated transfer function correction part 22 corrects the estimated transfer 

45 function by adding the estimated transfer function correcting vector 5h(k) to the estimated transfer function 
vector h(k) as expressed by the following equation. 

h(k+l) = b(k) + n8h(k) (4) 

50 

where u is called a step size, which is a preselected quantity for controlling the range of each correction 
and is handled as a time-invariant constant. In the following description, the estimated transfer function 
correcting vector ah(k) is calculated for it = 1 and, if necessary, it is multiplied by a desired step size u. In 
some applications the value of the step size ti is caused to vary with time, but in such a case, too, the 
55 following description is applicable. The corrected estimated transfer function vector is transferred to the 
convolution part 23. The above is a transfer function estimating operation at time k and the sam operation 
is repeated after time k + 1 as well. 
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The estimated transfer function correcting method described above in respect of Fig. 2 is known as an 
adaptive algorithm. While an LMS (Least Mean Square) algorithm and an NLMS (Normalized LMS) are also 
well-known as adaptive algorithms, a description will be given of a projection algorithm proposed in a 
literature [Ozeki and Umeda: "An Adaptive Filtering Algorithm Using an Orthogonal Projection to an Affine 

5 Subspace and Its Properties", Journal of Institute of Electronics, Information and Communication Engineers 
of Japan (A), J67-App. 126-132, (1984-2).] 

The projection algorithm requires a larger number of operations than does the NLMS but has an 
excellent adaptability to a speech signal input. With the projection algorithm, as referred to above, the 
vector 5h(k) is determined by Eq. (4) for u = 1 so that simultaneous equations composed of the following p 

to equations are satisfied. 

y(k) = x(k) T (h(k) + 5h(k) ) 
y(k-l) = x<k-l)T(h(k) + Sh(k)} 

75 

y(k-p+l) = x(k-p+l)T(h(k) + 5h(k)) (5) 

20 Eq. (5) indicates that the vector 5h(k) is determined so that the estimated transfer function h(k + 1) updated 
at time k provides the same values y(k), y(k-1), .... y(k + p-1) as the outputs from the unknown system, 
respectively, for p input vectors x(k), x(k-1), _ x(k + p-1) prior to time k. By this, it is expected that the 
characteristic of the estimated transfer function h(k) will approach the characteristic of the unknown system 
as the adaptive updating of the estimated transfer function is repeated using the vector 5h(k). 

25 In the above, p is a quantity commonly called a projection order. As the projection order p increases, 
the adaptability of the projection algorithm increases but the computational complexity also increases. The 
conventional NLMS method corresponds to the case of p = 1. 
Now, transposing the first equation in Eq. (5), we have 

30 

x(k> T Sh(k> = y(k) - x(k) T h(k) = e(k) (6) 

Furthermore, transposing the second equation in Eq. (5) and using an equation obtainable by setting k in 
35 Eqs. (4) and (6) to k-1 . we have 

x(k-l)T6h(k) = y(k-l) - x(k-l)Th(k) 
40 = y(k-l) - x(k-l)T ( h(k-l) + *iSh(k-l)) 

= y(k-l) - x(k-l)T ( h(k-i) + nx(k-l) T 5h(k-l) 
= e(k-l) - Jie(k-l) 

= (l-li)e(k-l) 



45 



50 



55 



(7) 

Thereafter, the following relationship similarly holds. 

x(k-i) T 5h(k) = (l-H)ie(k-i) (8) 

Based on this, Eq. (5) can be rewritten as the following system of simultaneous equations. 

X p (k)6h(k) = (k) (9) 

where Xp(k) is a matrix with p rows and L columns and e(k) is a vector of the p rows; they are defined by 
the following equations. 
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at (k)\ 
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(10) 
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(11) 



The vector e(k) will hereinafter be referred to as an error vector. Now, since p is usually smaller than L, Eq. 
to (9) is an under-determined simultaneous equation or indeterminate equation for the vector dh(k) and the 
minimum norm solution of the vector 5h(k) is given by the following equation. 



25 



5h(k) = X p (k) T (Xp(lc)Xp(k) T )- 1 e(k) 

= [x(k)x(k-D, x(k-lp+l) ]g(k) 



(12) 



where 

30 g(k) = R p (kpe(k) (13) 
Rp(k) = Xp(k)Xp(k) T (14) 

R p is a matrix with p rows and p columns, which will hereinafter referred to as a p-order covariance matrix 
35 or auto-correlation matrix, and g(k) a pre-filter coefficient vector. Letting elements of the pre-filter coefficient 
vector g(k) be represented by gi(k), g 2 (k), .... g p (k), the estimated transfer function correcting vector 5h(k) 
can be expressed on the basis of Eq. (12) as follows: 



40 



P 

5h (k) = I g. (k) x (k-i+1) 
i=l 1 



(15) 



45 Thus, when the projection algorithm is used, the estimated transfer function correcting vector calculation 
part 21 in Fig. 2 has such a construction as depicted in Fig. 3. Reference numeral 31 denotes a pre-filter 
coefficient vector calculation part, which uses the input signal x(k) and the error signal e(k) to calculate the 
pre-filter coefficient g(k) on the basis of Eq. (13). Reference numeral 32 denotes a pre-filtering part, which 
performs the pre-filtering operation expressed by Eq. (15) to synthesize the estimated transfer correcting 

so vector $h(k) by use of the pre-filter coefficient g(k) that is transferred from the pre-filter coefficient vector 
calculation part 31 . 

Now, a description will be given of the computational complexity of the projection scheme described 
above. The computational complexity mentioned herein is the number of multiplication-addition (or addition) 
operations necessary for estimating operations per unit discrete time. The computational complexity of Eq. 
55 (2) in the convolution part 23 of the tap number L in Fig. 2 is L. The computational complexity of Eq. (13) in 
the pre-filter coefficient vector calculation part 31 of the estimated transfer function correcting vector 
calculation part 21 is about p 3 ^ when using the Choleski method which is a typical computation method. 
The computational complexity of Eq. (15) in the pre-filtering part 32 is (p-1)L. The computational complexity 
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of Eq. (4) in the estimated transfer function correction part 22 is L. Thus, the entire computational 
complexity NC of the projection scheme is given by the following equation. 

NC = L + p 3 /6 + <p-1)L + L (16) 

5 

On the other hand, the computational complexity of the NLMS scheme or LMS algorithm is about 2L. 
For example, when L = 500 and p = 20 (a typical value in the case of an acoustic echo canceler), the 
number of operations involved in the NLMS scheme is 1000, whereas the projection scheme requires as 
many as about 12000 operations on the basis of Eq. (16). The computational complexity p 3 /6 of Eq. (13), in 
70 particular, abruptly increases as the projection order p increases. Thus, the projection scheme has excellent 
convergence characteristics as compared with the NLMS scheme but poses the problem of increased 
computational complexity. 

SUMMARY OF THE INVENTION 

15 

It is therefore an object of the present invention to provide an adaptive transfer function estimating 
method which permits reduction of the computational complexity involved and an estimating device using 
such a method. 

To attain the above objective, the present invention reduces the computational complexity of the pre- 
20 filter coefficient vector calculation part 31 according to its first aspect and the computational complexity of 
the pre-filtering part 32 according to its second aspect. 

According to the first aspect of the present invention, the adaptive transfer function estimating method 
or device which employs the projection algorithm and in which: 

the transfer function of an unknown system is estimated using an input signal x(k) and an output signal 
25 y(k) of the unknown system at a discrete time k; 

an error signal e(k) = y(k) - y(k) is calculated between an output signal y(k) of an estimated system 
having the estimated transfer function h(k) and.the output signal y(k) of the unknown system; 

letting the vector of the error signal e(k), a covariance matrix of the input signal x(k) and the vector of 
pre-filter coefficients be represented by e(k), R p (k) and g(k), respectively, the following simultaneous linear 
30 equation with p unknowns 

R P (k)g(k) = e(k) 

is solved to obtain the pre-filter coefficient vector g(k), and the pre-filter coefficient vector g(k) is used to 
35 calculate an estimated transfer function correcting vector 5h(k) by the following equation 

Sh(k) = [x(k) , x(k-l) x(k-l+p)]g(k); 

40 and the estimated transfer function correcting vector 5h(k) and a predetermined correcting step size u are 
used to repeatedly correct the estimated transfer function vector h(k) by the following equation 

h(k+l) = h(k) + Ji8h(k) 

45 

so that the error signal e(k) approaches zero; 

characterized in that a forward linear prediction coefficient vector a(k) of the input signal x(k), the sum 
of its a posteriori prediction-error squares F(k), a backward linear prediction coefficient vector b(k) of the 
input signal x(k) and the sum of its a posteriori prediction-error squares B(k) are computed and, letting a 
so pre-filter deriving coefficient vector be represented by f(k), the pre-filter coefficient vector g(k) is obtained 
by a recursive formula composed of the following first and second equations: 



55 
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According to the second aspect of the present invention, the adaptive transfer function estimating 
method or device which employs the projection algorithm and in which: 

the transfer function of an unknown system is estimated using an input signal x(k) and an output signal 
75 y(k) of the unknown system at a discrete time k; 

an error signal e(k) = y(k) - y(k) is calculated between an output signal y(kj of an estimated system 
having the estimated transfer function h(k) and the output signal y(k) of the unknown system; 

letting the vector of the error signal e(k), a covariance matrix of the input signal x(k) and the vector of a 
pre-filter coefficient be represented by'e(k), R p (k) and g(k), respectively, the following simultaneous linear 
20 equation with p unknowns 

R p (k)g(k) = e(k) 

is solved to obtain the pre-filter coefficient vector g(k). and the pre-filter coefficient vector g(k) is used to 
25 calculate an estimated transfer function correcting vector 5h(k) by the following equation 

Sh(k) - [x(k), x(k-l) x(k-l+p) ]g(k) ; 

30 and the estimated transfer function correcting vector 5h(k) and a predetermined correcting step size u are 
used to repeatedly correct the estimated transfer function vector h(k) by the following equation 

h(k+l) = h(k) + |i5h(k) 

so that the error signal e(k) approaches zero; 
characterized in that: 

instead of calculating said correcting vector 5h(k), pre-fifter coefficients g(k) which are elements of the 
pre-filter coefficient vector g(k) are smoothed by the following equation 

si(k) = si-i(k-l) .+ ^gi(k) for 2 < i < p 
= jigi (k) - for i = 1 



35 



40 



45 to obtain a smoothing coefficient Sj(k); 

instead of correcting said estimated transfer function vector, the smoothing coefficient s(k) is used to 
obtain an approximate estimated transfer function z(k + 1) by the following equation 

Z(k+1) = z(k) + s p (k)x(k-p + 1); 

50 

a convolution x(k) T z(k) between the approximate estimated transfer function z(k + 1) and the input signal 
x(k) is performed; 

an inner product Sp-i(k-l) T rp-i(k) is calculated, setting the vector of the smoothing coefficient s(k), the 
vector of the input signal x(k) and the correlation vector of the input signal as follows: 

55 

Sp-^k-l) = [si(k-1), s 2 (k-1), .... Sp-^k-l)] 7 

r P -i(k-1) = [x(k) T x(k-1), x(k) T x(k-2), .... x(k) T x(k-p + 1)] T ; and 
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the sum of the convolution result x(k) T z(k) and the inner product is output as the estimated signal y(k). 



BRIEF DESCRIPTION OF THE DRAWINGS 



5 Fig. 1 is a block diagram illustrating an ordinary construction for the estimation of a transfer function; 
Fig. 2 is a block diagram of a transfer function estimation part 11 in Fig. 1; 

Fig. 3 is a block diagram of an estimated transfer function correcting vector calculation part 21 
embodying the present invention in Fig. 2; 

Fig. 4 is a block diagram of a pre-filter coefficient vector calculation part 31 according to the first aspect 
10 of the invention in Fig. 3; 

Fig. 5 is a block diagram showing an estimated signal generating process according to the second 
aspect of the invention; 

Fig. 6 is a block diagram showing an application of the present invention to measurement of the transfer 
function of a loudspeaker; 
75 Fig. 7 is a block diagram showing an application of the present invention to an echo canceller; 
Fig. 8 is a block diagram showing an application of the present invention to noise control; 
Fig. 9 is a graph showing echo cancelling characteristics of an echo canceller embodying the present 
invention; and 

Fig. 10 is a graph showing, in comparison with a prior art example, the relationship between the number 
20 of operations for the estimation of the transfer function and the projection order. 

DESCRIPTION OF THE PREFERRED EMBODIMENTS 



A description will be given first of a method for reducing the computational complexity of the pre-filter 
25 coefficient vector calculation part 31. This method is characterized in that the pre-filter coefficient vector g- 
(k) is obtained by a recursion formula utilizing a vector g(k-1) at one unit time before. Now, it will be 
demonstrated how the vector g(k) can be obtained from the immediately preceding vector g(k-1) on a 
recursive basis. The recursion formula is derived through utilization of the fact that elements of the error 
vector e(k) of Eq. (11) shift with time while being multiplied by 1-u and the property of the inverse of the 
30 covariance matrix of the input signal given by Eq. (14). 

As is the case with Eq. (14), a covariance matrix R p -i(k) of the p-1 order is defined by the following 
equation. 



35 



40 



Vi (k) = Vi (kl Vi ,k) 



x(k) 
x (k-1)' 



x(k-p+2) J 



[x (k)x <k-l) (k-p+2) ] 

(17) 



45 



It is known from a literature, S. Haykin, "Adaptive Filter Theory," 2nd edition, Prentice-Hall, 1991, pp. 
577-578 that the inverse of the covariance matrix R p (k) and the inverse of the covariance matrix R p -i(k) 
bear such a relationship as shown below. 

50 



55 
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70 
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(19) 



where a(k) is a forward linear prediction coefficient vector (p-order) which satisfies a norma! equation R p (k)- 
a(k) = [F(k), 0, .... 0] T and its first element is a 1 ; F(k) is the minimum value of the sum of forward a 
posteriori prediction-error squares; b(k-1) is a backward linear prediction coefficient vector (p-order) which 

satisfies a normal equation R p (k-1)b(k-1) = [0 0, B(k-1)] T and its last element is a 1; and B(k-1) is the 

minimum value of the sum of backward a posteriori prediction-error squares. The linear prediction 
coefficients a(k), b(k-1) and the minimum values of the sums of forward and backward a posteriori 
prediction-error squares F(k) and B(k-1) could be calculated with less computational complexity through use 
of an FTF (Fast Transversal Filters) algorithm, for instance, as disclosed in J. M. Cioffi and T. Kailath, 
"Windowed fast transversal adaptive filter algorithms with normalization," IEEE Trans. Acoust, Speech 
Signal Processing, vol. ASSP-33, no. 3, pp. 607-625. The pre-filter coefficient vector g(k) of Eq. (13) is 
shown again below. 



g(k) = R P (k)- 1 e(k) (20) 

35 

A pre-filter deriving coefficient vector f(k) of p-1 order is defined, corresponding to the vector g(k), by the 
following equation. 

f(k) = R p -i(k)-Vi(k) (21) 

40 

where 



(22) 



50 



45 



e (k) = 

P-1 



e (k) 
(1- \i )e(k-l) 



(1- *i ) P ~ 2 e(k-p+2) 



Since all elements of a vector (1-u)e p _ l (k-1) of p-1 order, obtained by substituting (k-1) for k in Eq. (22) and 
multiplying both of its left and right sides by (1-u), constitute elements of the p-order vector (k) of Eq. 
(11) except its first element e(k), the relationship of the following equation holds. 



o 
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(k) = 



e(k) 



(l-H)e (H-l) 
P-1 



(23) 



Furthermore, since all the elements of the p-1 order vector of Eq. (22) constitute all elements of the p-order 
ro vector of Eq. (11) except its last element, the relationship between the vectors e(k) and ep-^k) can also be 
expressed as follows: 



75 



e (k) = 



(1- li ) 



. (k) 
p-1 

p-:L e(k-p+l) 



(24) 



20 



Multiplying the respective terms on both sides of Eq. (18) by e(k) from the right, we obtain the following 
equation from Eqs. (20), (21) and (23). 



25 



a (k) = (l- n ) 



|_f (k-l)J 



a (k) e (k) 
+ S — (k) (25) 



F(k) 



30 Moreover, multiplying both sides of Eq. (19) by e(k-1) from the right, we obtain the following equation from 
Eqs. (20), (21) and (24). 



«<*-!> .p"""!. «""-»'« <"-» b „c-i, ,26, 
|_ 0 j B(k-l) 



40 Transposing the right side to the left side, we have 



[ £ r]= 



T 

b (k-1) e (k-1) 

, n , ^(k-l) S — b(k-l) (27) 

45 L u B(k-l) 



In this way, the pre-filter coefficient vector g(k) is calculated on the basis of the value of the pre-filter 
deriving coefficient vector f(k-1) as expressed by Eq. (25), and the vector f(k-1) is calculated from the vector 

50 g(k-1) as expressed by Eq. (27). That is, Eqs. (25) and (27) are recursion formula for the vector g(k). 

On the right side of Eq. (25), to obtain the first term (1-u)[] needs p-1 multiplications and to obtain the 
second term a(k) T e(k) needs p multiply-add operations, one division for division by F(k) and p multiply-add 
operations for multiplication by a(k) and addition to the first term on the right side. That is, a total of 3p 
operations are needed. Similarly, approximately 2p operations are required to obtain Eq. (27). Accordingly, 

55 the number of operations of Eqs. (16) and (17) is around 5p, which is a computational complexity 
proportional to p. On the other hand, it is disclosed in the afore-mentioned literature by J. M. Cioffi and T. 
Kailath that the linear prediction coefficient vectors a(k), b(k-1) and the minimum values of the sum of a 
posteriori prediction-error squares F(k), B(k-1), which are needed in Eqs. (25) and (27), can be calculated by 
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w 



a linear prediction scheme with a computational complexity of about "lOp. Thus, it is evident that the 
recursion formulae of Eqs. (25) and (27) permits calculation of the vector g(k) with the computational 
complexity proportional to the projection order p. 

According to the pres nt invention based on the above-described principles, th matters (A), (B) and 
(C) listed below are effective in reducing the computational complexity and in increasing the computational 
stability. 

(A) Eq. (27) is an equation for the p-order vector and its left side is zero with respect to a p-th element 
(the least significant element) of the vector. Since a p-th element of the vector b{k-1) is always a 1, the 
following equation holds using a p-th element g p (k-1) of the vector g(k-1). 



b (k-l) T e 



(k-1) 



B(k-l) 



(28) 
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Therefore, the following equation may be calculated in place of Eq. (27). 
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[f (k-l)l 



1) -g (k-1) b (k-1) 
P 



(27a) 
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30 



In this case, however, taking computational errors into account, it may sometimes be advantageous to 
obtain both g p (k-1) and the left side of Eq. (28) and average them. 

(B) When denominators F(k) and B(k-1) of the second terms on the right sides of Eqs. (25) and (27) are 
small, operations are unstable. The computational instability could be reduced by adding non-negative 
5 F (k) and 5 B (k-1) to the denominators as shown in the following equations. 



35 



fir (k) = (i- 



f (k-l) 



(k) x e(k) 



F(k) 



5 p (k) 



(k) 



(29) 



40 




* (k-1) 



b (k-1) e (k-1) 

B(k-l) +5 (k-1) 
B 



b (k-1) 



(30) 



In practice, the values 5 F (k) and 5 B (k-1) may be set to desired values about 40 dB smaller than the 
45 average power of the input signal x(k) (that is, about 1/10000 the average power) or may also be varied 
with time k in accordance with power variations. 

(C) In linear prediction analysis for obtaining the linear prediction coefficient vectors a(k), b(k-1 ) and the 
minimum values of the sum of prediction-error squares F(k), B(k-1) which are needed in Eqs. (25) and 
(27), the analysis frame of the input signal x(k) differs from that in an ordinary case. In concrete terms, 

so the analysis frame in an ordinary linear prediction analysis ranges from time 0 to the current time, and 
when time k-1 is updated to k, x(k) needs only to be added to the analysis frame. In the projection 
algorithm, the analysis frame ranges from x(k) to x(k-L-p + 2) in Eq. (5); so that when time k-1 is updated 
to k, it is necessary not only to add x(k) to the analysis frame but also to remove x(k-L-p + 1). On this 
account, the linear prediction analysis in the projection algorithm requires computational complexity twice 

55 that needed for ordinary linear prediction analysis. However, when th statistical property of the input 
signal does not change or when it is expected that its change is slow, the result of the linear prediction 
analysis does not largely depend on the analysis frame; hence, it is possible to use the ordinary linear 
prediction analysis in which x(k) needs only to be added to the analysis frame at time k - this permits 
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reduction of the computational complexity. 

Next, a description will be given of an embodiment of the adaptive transfer function estimating method 
of the present invention based on the above-described theoretical discussions, reference being made to Fig. 
2 because the overall construction of its functional block is the same as shown in Fig. 2. Since the above- 

5 described theory underlying the present invention is related to the reduction of computational complexity in 
the estimated transfer function correcting vector calculation part 21 in Fig. 2, in particular, in the pre-filter 
coefficient vector calculation part 31 and the pre-filtering part 32 depicted in Fig. 3, these parts will 
hereinbelow be described in detail. 

Fig. 4 illustrates in block form the pre-filter coefficient vector calculation part 31 based on the 

w discussion made in the above. Reference numeral 41 denotes a linear prediction part, 42 a pre-filter 
deriving coefficient vector correcting part, 43 a pre-filter coefficient vector correcting part, and 44 an error 
vector generating part. The linear prediction part 41 calculates the forward linear prediction coefficient 

vector a(k) which satisfies a normal equation Rp(k)a(k) = [F(k), 0 0] T and the sum of forward a posteriori 

prediction-error squares F(k) when the prediction coefficient a(k), the backward linear prediction coefficient 

75 vector b(k) which satisfies a normal equation R p (k-1)b(k-1) = [0, .... 0, B(k-1)] T and the sum of backward a 
posteriori prediction-error squares B(k) when the prediction coefficient b(k) is used. These values can be 
calculated by methods disclosed in the afore-mentioned literature by J. M. Cioffi et al. The error signal 

vector generating part 44 stores p error signals e(k), e(k-1) t e(k-p-M) and constitutes the error signal 

vector e(k) of Eq. (11). In the pre-filter deriving coefficient vector correcting part 42, the pre-filter coefficient 

20 vector g(k-1), the backward linear prediction coefficient vector b(k-1), the error signal vector e(k-1) and the 
minimum value of the sum of backward a posteriori prediction-error squares B(k-1) are used to compute the 
pre-filter deriving coefficient vector f(k-1) on the basis of Eq. (27) or (30). In the pre-filter coefficient vector 
correcting part 43, the pre-filter deriving vector f(k-1), the error signal vector e(k), the forward linear 
prediction coefficient vector a(k) and the minimum value of the sum of forward a posteriori prediction-error 

25 squares F(k) are used to compute, by Eq. (25) or (29). the vector g(k) that satisfies Eq. (13). 

With the methods mentioned above, the computational complexity involved in the pre-filter coefficient 
vector calculation part 31 can substantially be reduced from p 3 /6 to 15p. Since the computational 
complexity (p-1)L in the pre-filtering part 32 remains unsolved, a large number of operations are still needed 
when L is large. 

30 Next, a description will be given of a solution to the above-noted problem by storing approximate values 
of the estimated transfer function vector h(k + 1) and averaging pre-filtering coefficients. 

At first, substituting k-1 for k in Eq. (4), the estimated transfer function h(k) is expressed as follows: 

35 h(k) = h(k-l) + llfih(k-l) (31) 

Substituting this in Eq. (4), the estimated transfer function h(k + 1) is expressed as follows: 
40 h(k+l) = h(k-l) + Jl5h(k) + ^Sh(k-l) (32) 

Setting k-2, k-3, ... for k in Eq. (4) and repeating the above substitution, the estimated transfer function 
h(k + 1 ) is expressed by 

h(k+l) = ji5h(k) + ^8h(k-l) + ... + |ifih(0) (33) 

In this instance, the estimated initial value fi(0) is set to 0 and this equation reveals that the estimated 

transfer function h(k + 1) is a summation of correcting vectors U5h(k), ufih(k-l) u6h(0) from time 0 (the 

so transfer function estimation starting time) to the current time k. 

The correcting vector u5h(k) is expressed by Eq. (15). Setting k-1, k-2, .... for k in Eq. (15) and 
substituting Eq. (15) in Eq. (33), we get 
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h(k+l) = p.[{gi(k)x(k)+g2(k)x(k-l) + ... + g p (k)x(k-p+l) } 
+ {gi(k-l)x(k-l)+g2(k-l)x(k-2} 
+ ... " +SfpfJc-l)x ? (k-p.) J 

♦ . . . • . r h * ; "• 

+ {gi (0)x(0)+g2(0)x(-l) + ... + g p (0)x(-p+l) } ] 
= jitgi(k)x(k) > : " 

+ '{g2(k>+gii(^ 

+ {g3(k)+g2(k^l)+g^(k-2^}x(k-2) , 1 

♦ • • ■ r ■ r, V • ' • ... 

+ {g p -l(k)+g p -2(k-l)+ . . . + gi(k-p+2) }x(k-p+2) 

+ (gp(k) +g p -i +. : . . + g2(k-p+2) 

+ gi(k-p+l)}x(k-p+l) + {g p (k-l)+g p -i(k-2) 

♦ ... + g2 (k-p+l)+gi(k-p) }x(k-p) 

♦ {gp(k-2)+g p -i(k-3)+ ... +g2(k-p) 

♦ gi<Jc-p-l)}x(k-p-l) 

♦ .. . (34) 

From this equation the following facts are expected. Firstly, the pre-filter coefficient gj(k) is calculated at 
every time k in the pro-Mtcf coefficient calculation part 31 in Fig. 3 and provided to the pre-filtering part 33; 
in this case, it »s expected that the computational complexity would be reduced by smoothing (or averaging) 
the pre-filter coefficient g,(k). Secondly, the term + g p (k-1)+ ... and the subsequent terms in Eq. (34) do not 
involve the pre-fitter coefficient gtfk) which is fed at time k, and hence this portion does not change after 
time k. By storing these terms as approximate values of the vector h(k + 1), it is expected that the 
computational complexity would be reduced accordingly, because no calculations are needed for these 
terms after time k 

Next, the above will be expressed by mathematical formulae. The smoothing of the pre-filter coefficient 
takes place for each corresponding input vector x(k-i). From Eq. (34), the smoothing corresponding to the 
input vector x(k-1). lor example, is g2(k) + gi(k-1), and the smoothing corresponding to the input vector x(k- 
2) is g3(k) + g:(k-t)*g»(k-2). Letting the result of averaging (a smoothing coefficient) corresponding to the 
input vector x(ki ♦ l). inclusive of the constant u, be represented by s { (k), it is expressed by 

i-1 

s.{k) = ii i g (k-j) far 1 < i < p (35) 
j=0 ^ 

i-1 

s.(k)= u y g. .(k-j) for p < i (36) 
i . . l-j 

D=i-p 



Eq. (35) is expressed as follows: 
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i-1 

s.(k) = n x g. . (k-j) + \ig. <k) 
j=l 1 

i-2 

= H X 9 i-j-i <lc ~ :i " 1) t ^i (lc) 



= s L _ 1 (k-l) ■■ + ^g i (k) - for 2 < i < p 

= ^g i (k) for i = 1 (37) 



On the other hand, letting the approximate value of the estimated transfer function h(k+1) be 
represented by z(k + 1), it is expressed by 

Z(k + 1) = {g P (k) + g p -,(k-1)+ ... 
+ 92<k-p + 2) + g, (k-p + 1 )}x(k-p + 1 ) 
+ {gp(k-1) + g P -,<k-2)+ ... 
+ g 2 (k-p + 1) + g,(k-p)}x(k-p) 
+ {9p(k-2)+g p -,(k-3)+ ... 
+ g2(k-p) + gi(k-p-1)}x(k-p-1) 
+ ... (38) 

From Eqs. (34), (35) and (38) we have the estimated transfer function fi(k + 1) expressed as follows: 



h (k+1) =z(k+l) + I s i (k)x (k-i+1) (39) 

i=l 



Moreover, the following relationship holds between the approximate values z(k + 1) and z(k). 



z(k+l) = {g p (k)+g p _i(k-l)+ . . . +g 2 (k-p+2) 
+gi(k-p+l) }x(k-p+l)+z(k) 
= s p (k)x(k-p+l) +z(k) (40) 

The estimated value y(k) of the output from the unknown system is expressed, from Eqs. (2) and 
(39), as follows: 
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y(k) = x(k) T h(k) 

P-l 

5 = x(k) T {z(k) + X si(k-l)x(k-i)} 

i=l 

P-l 

w = x(k)T z( k) + I .si(k-l)x(Jc)T x (k-i) 

i=l 

= x(k)T Z (k) + 8p-i(k-l)T rp -i(k) (41) 

75 

In the above, s p -i(k-1) is a smoothing coefficient vector and r^flO is a correlation vector, which are 
defined by the following equations. 

Sp-Uk-1) = [si(k-l), s 2 (k-1), Sp-^k-ljf (42) 

20 

r p -i(k) = [x(k) 7 x(k-1, x(k) T x(k-2). .... x(k) T x(k-p + 1 )] T (43) 
Since the vector x(k), defined by Eq. (3), is given by 

25 x(k) = [x(k). x(k-1 ) x(k-L + 1 )] T (44) 

the following relationship holds 

x(k) T x(k-i) = x(k-1 ) T x(k-i-1 ) -x(k-L)x(k-L-i) + x(k)x(k-i) (45) 

30 

where i = 1, 2, .... p-1 and the following equation holds 
r P -,(k) = rp-^k-l^k-DXp-^k-L) +x(k)x p . 1 (k) (46) 
35 where 

x P -,(k) = [x(k-1), x(k-2), .... x(k-p+ 1)] T (47) 

Next, a description will be given, with reference to Rg. 5, of the storage of the approximate value z(k) of 
40 the estimated transfer function h(k) and the transfer function estimation procedure which is followed when 
the pre-filter coefficients are smoothed. At time k, a correlation calculation part 52, which is supplied with 
the input signal x(k), calculates the correlation vector r p _!(k) by Eq. (46), using the input signal x(k), 

previous input values x(k-1) x(k-L) and the correlation vector r^k-l) at the immediately preceding 

time. 

45 Then, an inner product s p -i(k-1) T r p -i(k) of the correlation vector r p _ t (k) and the smoothing coefficient 
vector Sp-i(k-t) of the pre-filter coefficients is calculated in an inner product calculation part 53. A 
convolution part 54 performs a convolution x(k) T z(k) of the stored approximate value z(k) of the estimated 
transfer function and the input signal. The results of the inner product calculation and the convolution are 
added together in an addition part 57 to synthesize the estimated value y(k) of the unknown system output. 

so These operations correspond to the operation of Eq. (41). 

Following this, the error e(k) between the unknown system output y(k) and the estimated value y(k) is 
obtained by the subtractor 24 shown in Rg. 2 and the pre-filter coefficient gj(k) is calculated in the pre-filter 
coefficient vector calculation part 31 shown in Rg. 4. 

After this, the calculated pre-filter coefficients are sent to a pre-filter coefficient smoothing part 51 in 

55 Rg. 5, wherein they are smoothed to obtain p smoothing coefficients Si(k), S2(k), .... Sp-i(k), s p (k). This 
smoothing operation is performed on the basis of Eq. (37). Of the smoothing coefficients, Si(k), S2(k), .... 
Sp-^k) are supplied, as elements of th smoothing coefficient vector Sp-i(k), to the inner product calculation 
part 53 and an estimated transfer function calculation part 56. Th smoothing coefficient s p (k) is fed to an 
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estimated transfer function approximate value storage part 55. 

The estimated transfer function approximate value storage part 55 updates the approximate value, using 
the smoothing coefficient s p (k) and the input signal vector x(k). That is, s p (k)x(k-p + 1 ) is added to the 
approximate value z(k) stored until then and the added result is stored as z(k + 1). This operation 
5 corresponds to the computation Eq. (40). 

Finally, in the estimated transfer calculation part 56 the smoothing coefficient vector Sp-^k) and the 
input signal vector x(k) are used to calculate Eq. (39) to obtain the estimated transfer function h(k + 1). 

In the above-described operations, rough estimates of computational complexity involved in the 
respective parts are as follows: 
10 Correlation calculation part 52: 2p , 

Inner product calculation part 53: p 

Convolution part 54: L 

Pre-filter coefficient vector calculation part 31 : 1 5p 
Pre-filter coefficient smoothing part 51: p 
75 Linear approximate value storage part 55: L 

Estimated transfer function calculation part 56: Lp 
In the above, p-1 is regarded as nearly equal to p. The entire computational complexity is such as follows: 

(L+19)p + 2L (48) 

20 

Now, attention should be paid to the following points. With the conventional transfer function estimation 
method depicted in Fig. 2, the estimated transfer function h(k + 1) corresponding to the transfer function is 
calculated at every time and is used to synthesize the estimated value y(k) of the unknown system output. 
In contrast thereto, according to the present embodiment which utilizes the approximate value z(k) of the 

25 estimated transfer function, the estimated value y(k) of the unknown system output can be synthesized 
without the need of calculating, directly, the estimated transfer function h(k + 1) as shown in Fig. 5. Once the 
estimated value y(k) is obtained, the above-described operations can be performed. Accordingly, there is no - 
need of calculating the estimated transfer function at every time in the cases where the estimated value 
h(k + l) of the transfer function, obtained as the result of the long-time estimating operation, is needed and 

30 where the estimated value y(k) of the unknown system output is needed rather than the estimated result of 
the transfer function (for example, in the case of the estimation of characteristics of a time-invarient system 
or in an acoustic echo canceller). 

On this account, tho overall computational complexity at every time, except the computational complex- 
ity in the estimated transfer function calculation part 56, is as follows: 

35 

19p + 2L (49) 

As referred to previously with respect to the prior art, the number of operations needed in the past is about 
12000 when the lap number L of the filter is 500 and the projection order p is 20. From Eq. (49), however, 

40 the number of operations m the present invention is 1380; that is, the computational complexity is reduced 
down to about 1/8 that m the prior art. 

As described above, the present invention permits substantial reduction of the computational complexity 
involved in the conventional projection scheme by the recursive synthesis of the pre-filter coefficient, 
storage of approximate values of the estimated transfer function vector and smoothing of the pre-filtering 

45 coefficients. 

While in the above the step size u has been described to be handled as a scalar, it may also be 
provided as uA by use of a diagonal matrix A. For example, in the case where the energy of an impulse 
response of an unknown system decays exponentially, the transfer function estimation speed may 
sometimes be increased by arranging the elements of the diagonal matrix A as a progression which decays 
so exponentially, as shown below. 

A = diag(a, a 7 , oy 2 a 7 L ~ 1 ), (0 < 7 < 1) (50) 

In this instance, the input to the linear prediction part is multiplied by a ratio 7 1/2 and R p (k) Eq. (14) is 
55 redefined by the following equation. 
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R p (k) = 



X (k) 1 
x (k-l) T: ' 



x (k-p+1) 



A tx(k)x(k-l) ... x (k+p-1) (51) 



io This correction permits the application of the present invention described above. 

Next, a description will be given of examples of application of the transfer function estimating device of 
the present invention. 

A first example of application is the measurements of transfer functions of pieces of acoustic 
equipment. Fig. 6 illustrates a system for measuring the transfer function from a loudspeaker to a 

is microphone. In Fig. 6, reference numeral 121 denotes a loudspeaker, 122 a microphone and 11 a transfer 
function estimating device. The output signal y(k) of the microphone 122 is a signal that has the 
characteristics of the loudspeaker 121 added to the input signal 1 x(k). Regarding the loudspeaker 121 
(including an acoustic path and the microphone as well) as an unknown system, the illustrated system is the 
same as that shown in Fig. 1, and by connecting the transfer function estimating device 11 of the present 

20 invention to the input of the loudspeaker 121 and the output of the microphone 122 as shown in Fig. 6, the 
transfer function h(k) of the loudspeaker can be estimated as the filter coefficient h(k) of the FIR filter. The 
measurement is usually performed in an anechoic chamber to avoid the influence of room acoustical 
characteristics. 

A second example is an acoustic echo canceller which prevents howling or echo in a loudspeaking 

25 communication system such as a TV conference system or visual telephone. Fig. 7 is a diagram for 
explaining the acoustic echo canceller. In Fig. 7, reference numeral 121 denotes a receiving loudspeaker, 
122 a transmitting microphone, 123 a room acoustic system and 20 an estimated echo; generating part 
which is identical in construction with the estimated signal generating part 20 embodying the present 
invention in Fig. 2. The transfer function estimating device 1 1 of this embodiment operates as an acoustic 

30 echo canceller. In a hands-free communication system using the loudspeaker and the microphone, the other 
party's voice emanating from the receiving loudspeaker 121 is received by the transmitting loudspeaker via 
the room acoustic system of the transfer function h(k). The received signal y(k) is returned to the other 
party and reproduced there. At the other party side, the transmitted voice is sent back and reproduced; this 
phenomenon is called an acoustic echo and disturbs comfortable communication. 

35 The estimated echo generating part 20 of the acoustic echo canceller 1 1 estimates the room acoustical 
characteristics h(k) including characteristics of the loudspeaker and generates an estimated echo signal y- 
(k) based on the characteristics h(k) estimated as those of the input signal x(k). The subtractor 24 subtracts 
the estimated echo signal y(k) from the received signal y(k). When the estimation is performed well, the 
echo canceller 1 1 operates so that it minimizes the power of the error signal e(k) and hence makes the 

40 estimated echo signal y(k) nearly equal to the received signal y(k), substantially reducing the acoustic echo. 
Comparison of the Fig. 7 system with the Fig. 2 system reveals that the transfer function estimating 
device of the present invention can be used directly as an echo canceller. The room acoustic system 123 in 
Fig. 7 corresponds to the unknown system 12 in Fig. 2 and the estimated echo y(k) and the transmitted 
signal e(k) in Fig. 7 correspond to the output y(k) from the convolution part 23 and the error signal e(k) in 

45 Fig. 2, respectively. 

A third example is noise control. Fig. 8 illustrates the principles of noise control. In Fig. 8, reference 
numeral 131 denotes a noise source, 132 a noise transfer path expressed by the transfer function h(k), 133 
a microphone for observation, 134 a noise monitor microphone, 20 an estimated noise generating part, 136 
a phase inverter and 137 a loudspeaker. The purpose of noise control is to cancel noise that is observed via 

so the noise transfer path of the transfer characteristics h(k) from the noise source 131, by g nerating negative 
estimated noise (letting a sound pressure be represented by y(k), a sound pressure represented by -y(k) is 
called a negative sound with respect to y(k)) from the loudspeaker 137. 

To attain this object, the transfer characteristics h(k) of the noise transfer path 132 is estimated in the 
estimated noise generating part 20 of the same construction as that of the estimated signal generating part 

55 20 in Fig. 2. That is, noise is detected by the monitor microphone 134 in the vicinity of the noise source 
131 and provided as the input signal x(k) to the estimated noise generating part 20, which generates an 
estimated value y(k) of the noise signal y(k) at the obs rvation point (i.e. the microphone 133). The 
estimated noise y(k) is polarity inverted by the phase inverter 136 into a signal -y(k). Assuming, for the sake 
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of brevity, that the loudspeaker characteristics are negligible, the signal -y(k) from the loudspeaker 137 is a 
combination of the noise signal y(k) and the sound pressure in the microphone 133 at the observation point 
and the error signal e(k) is provided at the output of the microphone 133. In this instance, when the noise 
transfer path is estimated well, the signal y(k) becomes similar to the noise y(k) from the noise source 131 
5 and the noise y(k) is cancelled by the combined sound pressure -y(k) from the loudspeaker 137. 

The microphone 133 in Fig. 8 corresponds to the subtractor 24 in Fig. 2; accordingly, the systems of 
Figs. 8 and 2 differ only in whether the error signal generated outside the estimating device 1 1 is provided 
thereto or the error signal is calculated in the estimating device 11. The principles of the present invention 
are applicable to the noise control device of Fig. 8. 

70 In the measurement of the transfer characteristics of the loudspeaker shown in Fig. 6, the speaker 
characteristics do not change with time; hence, there is no need of learning measured results (transfer 
characteristics) in the course of measurement and the transfer characteristics need only to be made known 
as the result of measurement conducted after a certain period of time. In the acoustic echo canceller shown 
in Fig. 7, the room transfer function varies with time in response to movement of audiences or opening and 

75 closing of a door, for instance. As is evident from Fig. 7, however, the purpose of the acoustic echo 
canceller is attained by obtaining the estimated value y(k) of the output y(k) of the room transfer function 
(the unknown system). Therefore, the estimated value h(k) of the room transfer function itself is unnec- 
essary in the acoustic echo cancejler. Also in the noise control device of Fig. 8, its purpose can be attained 
by obtaining the estimated value y(k) of the unknown system and no estimated value of the room transfer 

20 function itself is needed. 

As will be appreciated from the above, the transfer function estimating device of the present invention, 
when applied as described above, permits substantial reduction of the computational complexity as 
compared with the prior art. Furthermore, the present invention has a basic feature of minimizing the power 
of the error signal e(k) in such a system configuration as shown in Fig. 2. Thus, the present invention is 

25 applicable to all instances in which problems to be solved can be modelled as the error signal power 
minimizing problem shown in Fig. 2. 

In the embodiments described above, it will be effective to select the elements of the diagonal matrix A 
of Eq. (50) so that they provide the same attenuation as the room reverberation. 

Finally, a description will be given of experimental results of the acoustic echo canceller according to 

30 the present invention. Fig. 9 shows what are called learning curves, the ordinate representing echo return 
loss enhancement (hereinafter referred to as ERLE) and the abscissa time. As the estimation of the room 
transfer function proceeds with the lapse of time, the ERLE increases. In Fig. 9, the curve 21 1 is a learning 
curve in the case of the projection order p being 2, the curve 212 a learning curve in the case of the 
projection order p being 8 and the curve 213 a learning curve in the case of the projection order p being 32. 

35 It will be understood, from Fig. 9, that the larger the projection order, the higher the convergence speed (the 
ERLE increases in a short period of time). 

Fig. 10 is a graph showing the relationship between the projection order p and the required computa- 
tional complexity, the abscissa representing the projection order p and the ordinate the number of multiply- 
add operations (including add operations as well). In Fig. 10, the curve 221 indicates the computational 

40 complexity in the prior art and the curve 222 the computational complexity when the present invention was 
used. The tap number L of the FIR filter for use in the estimation of the transfer function was 500. From Fig. 
10 it will be seen that the computational complexity involved in the present invention is particularly smaller 
than in the case of the prior art when the projection order p is large. 

As described above, the present invention permits substantial reduction of the number of operations 

45 necessary for the estimation of the transfer function or output of an unknown system by the projection 
scheme. In concrete terms, letting the tap length of the FIR filter which represents the unknown system and 
the projection order be represented by L and p, respectively, the prior art requires multiply-add operations 
p 3 /6 + (p + 1)L, that is, the computational complexity is proportional to p 3 , whereas according to the present 
invention the computation complexity can be reduced down to I9p + 2L. 

so Such a reduction of computational complexity allows corresponding reduction in the scale of hardware, 
and hence significantly contributes to the downsizing of the device and reduction of its cost. Moreover, 
when the hardware is held on the same scale, the tap length L and the projection order can be chosen 
larger than in the past - this speeds up the estimating operation and increases the estimation accuracy. 
Besides, when the present invention is embodiment by a computer, the operation time can be greatly 

55 shortened. 

It will be apparent that many modifications and variations may be effected without departing from the 
scope of the novel concepts of the present invention. 
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Claims 

1. A method for adaptive estimation of a transfer function by use of a projection algorithm, said method 
comprising: .i ■ 

a step wherein the transfer functovof an -unknown system is estimated using an input signal x(k) 
and an output signal y(k) of said unknown 1 system at a discrete time k; ' * 

a step wherein an error signal &(ky "= ' y(k)^y(k) between an output signal y(k) of an estimated 
system having said estimated transfer function h(k) and said output signal y(k) of said unknown system 
is calculated; 

a step wherein, letting the vector of said error signale(k), a covariance matrix of said input signal x- 
(k) and the vector of a pre-filter coefficient be represented by e(k); Rp(k) and g(k), respectively, the 
following simultaneous linear equation with p unknowns ; 

R P (k)g(k) = e(k) 

is solved to obtain said pre-firter coefficient vector g(k), and, letting the vector of said input signal x(k) 
and an estimated transfer function correcting vector be represented by x(k) and 5h(k), respectively, said 
pre-firter coefficient vector g(k) is used to calculate said correcting vector 5h(k) for said estimated 
transfer function of said estimated system by the following equation 

6h(k) = [x(k), x(k-l), . .., x(k-l+p) ]g(k) ; 

and 

a step wherein said estimated transfer function correcting vector 5h(k) and a predetermined 
correcting step size u are used to repeatedly correct said estimated transfer function vector h(k) by the 
following equation - 

h(k+l) = b(k) + ji&Mk) 

so that said error signal e(k) approaches zero; 

wherein said step of obtaining said pre-firter coefficient vector g(k) by solving said simultaneous 
linear equation with p unknowns is a step wherein (a) a forward linear prediction coefficient vector a(k) 
of said input signal x(k), the sum of its a posteriori prediction-error squares F(k), a backward linear 
prediction coefficient vector b(k) of said input signal x(k) and the sum of its a posteriori prediction-error 
squares B(k) are computed and, (b), letting a pre-filter deriving coefficient vector be represented by f- 
(k), said pre-filter coefficient vector g(k) is obtained by a recursion formula composed of the following 
first and second equations: 

[0 ~| a (k) T e(k) 
— 7JZ\ — a(k) 
f(k-l)J F(k) 

P'^'l-ack-i) - b ' k - 1|T " k - 1 ' b , k -i l 
L 0 J B(k-l) 

2. A method for adaptive estimation of a transfer function by use of a projection algorithm, said method 
comprising steps: 

the transfer function of an unknown system is estimated using an input signal x(k) and an output 
signal y(k) of said unknown system at a discrete time k; 

an rror signal e(k) = y(k) - y(k) between an output signal y(k) of an estimated system having said 
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estimated transfer function h(k) and said output signal y(k) of the unknown system is calculated; 

letting the vector of said error signal e(k), a covariance matrix of said input signal x(k) and the 
vector of a pre-filter coefficient be represented by e(k), R p (k) and g(k), respectively, the following 
simultaneous linear equation with p unknowns 

R P (k)g(k) = e(k) 

is solved to obtain said pre-filter coefficient vector g(k), and said pre-filter coefficient vector g(k) is used 
to calculate a correcting vector 5h(k) for said estimated transfer function of said estimated system by 
the following equation 

5h(k) = [x(k) , x(k-l), . . . , x(k-l+p) }g(k) ; 

and said estimated transfer function correcting vector 5h(k) and a predetermined correcting step size u 
are used to repeatedly correct said estimated transfer function vector h(k) of said estimated system by 
the following equation 

h(k+l) = h(k) + Ji5h(k) 

so that said error signal e(k) approaches zero; 
characterized by the steps: 

instead of calculating said correcting vector 5h(k), pre-filter coefficients g(k) which are elements of 
said pre-filter coefficient vector g(k), are smoothed by the following equation 

si(k) = si-i(k-l) + \igi(k) for 2 < i < p 
= ^gi(k) for i = 1 

to obtain a smoothing coefficient Sj(k); 

instead of correcting said estimated transfer function vector by said correcting vector, said 
smoothing coefficient s(k) is used to obtain an approximate estimated transfer function z(k + 1) by the 
following equation 

z(k + 1) = z(k) + s p (k)x(k-p + 1); 

a convolution x(k) T z(k) of said approximate estimated transfer function z(k + 1) and said input signal 
x(k) is performed; 

wherein an inner product Sp-^k-l^rp-^k) is calculated, setting the vector of said smoothing 
coefficient s(k) and the correlation vector of said input signal to 

Sp-i(k-l) = [s,(k-1),s 2 (k-1) s p -,(k-1)] T 

and 

r P -i(k) = [x(k) T x(k-1) T x(k-2), x(k) T x(k-p + 1)] T ; and 

wherein the sum of said convolution result x(k) T z(k) and said inner product is output as said 
estimated signal y(k). 

The method of claim 1, further comprising a step wherein pre-filter coefficients g(k), which are elements 
of said pre-filter coefficient vector g(k), are smoothed by the following equation 

si(k) = si-i<k-l) + Hgi(k) for 2 < i < p 
= ^gi(k) for i = 1 

on 
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to obtain a smoothing coefficient Sj(k); 

a step wherein said smoothing coefficient s(k) is used to obtain an approximate estimated transfer 
function z(k + 1) by the following equation 

s z(k + 1) = z(k) + s p (k)x<k-p + 1); 

a step wherein a convolution x(k) T z(k) of said approximate estimated transfer function z(k + 1) and 
said input signal x(k); 

a step wherein an inner product s p -i(k-1) T r p -i(k) is calculated, setting the vector of said smoothing 
w coefficient s(k) and the correlation vector of said input signal to 

s p -i(k-1) = [si(k-1), s 2 (k-1), .... Sp-^k-l)] 7 

and 

rp-i(k) = [x(k) T x(k-1), x(k) T x(k-2) x(k) T x(k-p + 1)] T ; and 

a step wherein the sum of said convolution result x(k) T z(k) and said inner product is output as said 
estimated signal y(k). 

4. The method of claim 1 or 3, wherein, letting the last element of said pre-filter coefficient vector g(k-t) 
be represented by g p (k-1), said pre-filter coefficient vector g(k) is calculated using the following 
equation which is a modified version of said second equation: 
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[f (lc-l)l 



1) - g (k-l) b (k-l) 
P 



5. The method of claim 1 or 3, wherein, letting predetermined two non-negative numbers be represented 
by 5 F (k) and i B (k-l). said pre-filter coefficient vector g(k) is calculated by the following equations which 
are modified versions of said first and second equations: 



[0 ~l a (k) T e(1 
f(k-l>J + F(k> + 5 F 



*(k) = (i-n)| ; h ^rrr^ »< k > 



[ f r]= 



- b <■»-»*. <k-l) 

* 0 B(Jc-l) + 5 (k-l) 

" J B 



so 6. A device for adaptive estimation of a transfer function by a projection algorithm, which comprises: 

an estimated system formed by convolution means which has a transfer function h(k) estimated 
from an input signal x(k) to and an output y(k) from an unknown system at a discrete time k; 

pre-filter coefficient vector calculating means whereby, letting the vector of an error signal e(k) = 
y(k) -y(k) between an output signal y(k) of said estimated system and said output y(k) of said unknown 
55 system, a covariance matrix of said input signal x(k) and a pre-filter coefficient vector be represented 
by e(k), Rp(k) and g(k), respectively, the following simultaneous linear equation with p unknowns is 
solved to obtain said pre-filter coefficient vector g(k), 
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Rp(k)9(k) = (k); 

pre-filtering means which, letting the vector of said input signal x(k) and an estimated transfer 
function correcting vector be represesnted by x(k) and 5h(k), respectively, uses said pre-filter coefficient 
vector g(k) to calculate said correcting vector 5h(k) for said estimated transfer function of said estimated 
system by 

fih(k) = [x(k), x(k-l), . x(k-p+l) ]g(k) ; 

and 

estimated transfer function correcting means for correcting said estimated transfer function vector 
h(k) of said estimated system through use of said estimated transfer function correcting vector 5h(k) 
and a predetermined correcting step size u by 

h(k+l) = h(k) + ^Sh(k) ; 

and 

in which said estimated transfer function vector h(k) is repeatedly corrected so that said error signal 
e(k) approaches zero; 

wherein said pre-filter coefficient vector calculating means comprises: 

linear prediction means for calculating a forward linear prediction coefficient vector a(k) of said 
input signal x(k), the sum of forward a posteriori prediction-error squares F(k), a backward linear 
prediction coefficient vector b(k) and the sum of backward a posteriori prediction-error squares B(k); 

pre-filter coefficient vector correcting means whereby, letting a pre-filter deriving coefficient vector 
f(k), said pre-filter coefficient vector g(k) is corrected by a first one of the following first and second 
equations which forms a recursion formula for said pre-filter coefficient vector, 

[0 ! a (k) T e( k) 
^ — a<ki 



pr]- 



* (k -l) - b(k-l) T «(ic-i) 
B(k-l) 



and 

pre-filter deriving coefficient vector correcting means for correcting said pre-filter deriving coeffi- 
cient vector f(k) by said second equation. 

A device for adaptive estimation of a transfer function by a projection algorithm, which comprises: 

an estimated system formed by convolution means which has a transfer function h(k) estimated 
from an input signal x(k) to and an output y(k) from an unknown system at a discrete time k; 

pre-filter coefficient vector calculating means whereby, letting the vector of an error signal e(k) = 
y(k)-y(k) between an output signal y(k) of said estimated system and said output y(k) of said unknown 
system, a covariance matrix of said input signal x(k) and a pre-filter coefficient vector be represented 
by e(k), R p (k) and g(k), respectively, the following simultaneous linear equation with p unknowns is 
solved to obtain said pre-filter coefficient vector g(k) f 

Rp(k)g(k) = (k); 

pre-filtering means which, letting the vector of said input signal x(k) and an estimated transfer 
function correcting vector be represented by x(k) and 5h(k), respectively, uses pre-filter coefficient 

oo 
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vector g(k) to calculate said correcting vector 6h(k) for said estimated transfer function of said estimated 
system by 

5h(k) = (x(k), x(k-l) x(k-p+l) ]g(k) ; 

and estimated transfer function correcting means for correcting said estimated transfer function 
vector h(k) of said estimated system through use of said estimated transfer function correcting vector 
5h(k) and a predetermined correcting step u by 

70 

h(k)<k+l) = h(k) + llSh(k); 

and 

in which said estimated transfer function h(k) is repeatedly corrected so that said error signal e(k) 
approaches zero; 

wherein instead of calculating said correcting vector 6h(k), said pre-filtering means is pre-filter 
coefficient smoothing means for smoothing pre-filter coefficients g(k), which are elements of said pre- 
filter coefficient vector g(k), by 

si(k) = si-i(k-l) + Hgi(k) for 2 < i < p 
= ^gi(k) for i = 1 

and for providing said smoothed output as a smoothing coefficient srfk); and 
wherein said estimated transfer function correcting means comprises: 

approximate estimated transfer function calculating means which, instead, of correcting said * ■ = 
estimated transfer function vector by said correcting vector, uses said smoothing coefficient s(k) to 
obtain an approximate estimated transfer function z(k + 1) by 

z(k + 1) = z(k) + s p (k)x(k-p + 1); 

convolution means which performs a convolution x(k) T z(k) of said approximate estimated transfer 
function z(k + 1) and said input signal x(k); 
35 inner product calculating means which calculates an inner product 

s p -i(k-1)Vi(k), 

setting the vector of said smoothing coefficient s(k) and the correlation vector of said input signal x(k) to 

40 

Sp_i(k-1) = [si(k-1), s 2 (k-1) Sp-nfk-l)] 1 

and 

45 fp-^k) = [x(k) T x(k-1), x(k) T x(k-2) x(k) T x(k-p + 1)f; and 

adding means for providing, as said estimated signal y(k), the sum of the result of said convolution 
x(k) T z(k) and said calculated inner product. 

50 8. The device of claim 7, wherein said inner product calculating means includes correlation calculating 
means which, letting the tap number of said convolution means be represented by L, calculates said 
correlation vector of said input signal x(k) by 

rp-i(k) = rp-^k-lhxCk-LJvUk-g + tkJviW 

55 

where 

Xp-^) = [x(k-1), x(k-2) x(k-p + 1)] T 

9^ 
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which is the vector of said input signal x(k). 

9. The device of claim 6, wherein said pr -filtering means is pre-filter coefficient smoothing means for 
smoothing pre-filter coefficients g(k), which are elements of said pre-filter coefficient vector g(k), by 



and for providing said smoothed output as a smoothing coefficient s 4 (k), and wherein said estimated 
transfer function correcting means comprises: 

approximate estimated transfer function calculating means which uses said smoothing coefficient s- 
(k) to obtain an approximate estimated transfer function z(k + 1) by 

z(k + 1) = z(k) + s p (k)x(k-p + 1); 

convolution means which performs a convolution x(k) T z(k) of said approximate estimated transfer 
function z(k + 1) and said input signal x(k); 

inner product calculating means which calculates an inner product 

Sp-,(k-1)Vi(k). 

setting the vector of said smoothing coefficient s(k) and the correlation vector of said input signal x(k) to 
s p -i(k-1) = [Mk-1), s 2 (k-1), .... s p ^(k-1)] T 



rp-^k) = [x(k) T x(k-1), x(k) T x(k-2) x(k) T x(k-p + 1)] T ; and 

adding means for providing, as said estimated signal y(k), the sum of the result of said convolution 
x(k) T z(k) and said calculated inner product. 

10. The device of claim 6 or 9, wherein said pre-filter deriving coefficient vector correcting means is means 
which, letting the last element of said pre-filter coefficient vector g(k-1) be represented by g p (k-1)- 
,calculates said pre-filter coefficient vector g(k) by the following equation which is a modified version of 
said second equation: 



11. The device of claim 6 or 9, wherein, letting predetermined non-negative numbers be represented by fir 
(k) and 5 B (k-1), respectively, said pre-filter coefficient vector correcting means and said pre-filter 
deriving coefficient vector correcting means calculate said pre-filter coefficient vector g(k) by the 
following equations, respectively, which are modified versions of said first and second equations: 



si(k) = Si-i(k-l) + Hgi(lc) 



for 2 < i < p 
for i = 1 



and 
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[0 ~1 a (k) T e 
f(k-l)J FCW + 5 



*<*> = u-iDi - i* ,-: - - 8 ( F k (L) a(k) 



[ f r]= 



_ . (t .„. brt-u'.tfc-l) 

10 10 1 B(k-l) + 5 (k-1) 

B 



75 12. The device of claim 6 or 9, further comprising subtracting means for providing, as said error signal e(k), 
the difference y(k)-y(k) between said output signal y(k) of said unknown system and said estimated 
signal y(k) which is the output of said estimated system. 
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